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ABSTRACT 

We present a Monte-Carlo model of He II reionization by QSOs and its effect on the 
thermal state of the clumpy intergalactic medium (IGM). The model assumes that 
patchy reionization develops as a result of the discrete distribution of QSOs. It includes 
various recipes for the propagation of the ionizing photons, and treats photo-heating 
self-consistently. The model predicts the fraction of He III, the mean temperature in the 
IGM, and the He II mean optical depth — all as a function of redshift. It also predicts 
the evolution of the local temperature versus density relation during reionization. Our 
findings are as follows: The fraction of He III increases gradually until it becomes close 
to unity at z ~ 2.8 — 3.0. The He II mean optical depth decreases from r ^ 10 at 
z > 3.5 to T < 0.5 at z < 2.5. The mean temperature rises gradually between z ~ 4 
and z ~ 3 and declines slowly at lower redshifts. The model predicts a flattening of the 
temperature-density relation with significant increase in the scatter during reionization 
at z ~ 3. Towards the end of reionization the scatter is reduced and a tight relation 
is re-established. This scatter should be incorporated in the analysis of the Lya forest 
at z < 3. Comparison with observational results of the optical depth and the mean 
temperature at moderate redshifts constrains several key physical parameters. 

Key words: intergalactic medium - cosmology: theory - quasars - dark matter - large 
scale structure of the universe 



1 INTRODUCTION 

The physical properties and the thermal history of the 
intergalactic- medium (IGM), play an important role in shap- 
ing the observed structure in the universe. A proper mod- 
elling of the evolution of the IGM is therefore essential for 
the interpretation of the rapidly accumulating data on the 
high redshift universe. These data include the secondary 
temperature and polarization fluctuations in high resolution 
measurements of the cosmic microwave background (CMB), 
and absorption features in QSO spectra, which probe the 
ionized and neutral gaseous components, respectively. Since 
the distribution of gas follows the underlying mass, these 
data can be used to constrain the linear mass power spec- 
trum (Croft; et al. 1998, 2002; McDonald et al. 2000, 2004; 
Viel, Haehnelt & Springel 2004). 

The gaseous component of the IGM follows, by large, 
the distribution of the underlying mass distribution of the 
gravitationally dominant dark matter (DM). The temper- 
ature, ionization state and other physical properties are. 



however, determined by energetic feedback incurred by QSO 
and galactic activities. This feedback is manifested in photo- 
ionization and heating of the IGM by radiation emitted by 
QSOs and galaxies, and in mechanical energy from super- 
novae explosions of massive stars produce galactic winds 
that can shock-heat the IGM and enrich it with metals (e.g. 
Cowie & Songaila 1998; Aracil et al. 2004). 

Absorption features in QSO spectra indicate that me- 
chanical energy in galactic winds stirs up the IGM only 
in the vicinity of galaxies. Observations (Adelberger et al. 
2003) have shown that winds can evacuate the neutral hy- 
drogen from within a distance of up to 1-2 h~^Mpc away 
from Lyman break galaxies (LBGs). Away from galaxies, 
analysis of the Lya forest strongly suggests that the mod- 
erate density IGM which makes up most of the volume in 
space, is only affected by ionizing UV photons emitted by 
QSOs and galaxies. 

The gaseous IGM is mainly made of a primordial com- 
position of hydrogen and helium with the former being the 
dominant element. At z < 1000 these elements recombined 
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and remained neutral until the appearance of stars (and pos- 
sibly very high redshift QSOs) able to produce sufficient 
photons to cause significant reionization. Measurement of 
the polarization anisotropies of the CMB by the Wilkinson 
Microwave Anisotropy Probe (WMAP) satellite call for an 
ionized IGM at redshifts as high as 2; ~ 17 ± 5 (Kogut et 
al. 2003; Spergel et al. 2003). The Gunn-Peterson trough in 
QSO spectra seems to occur at 2 « 6 (Becker et al. 2001; 
Djorgovski et al. 2001) and might imply a complex reion- 
ization history of the IGM with more than a single phase 
of reionization (Gen 2003; Wyithe & Loeb 2003b). Ionizing 
radiation from galaxies is also sufficient to singly ionize he- 
lium (Benson et al. 2002) at 2: > 6. In any case observations 
of the Lya forest show that by z ~ 6 most of the hydro- 
gen in the universe was photo-ionized and that the IGM is 
in a quiet state governed by the gravitational drag of the 
dark matter and the photo-heating of the ionizing radiation 
(see Rauch 1998 for a review). During this stage a tight re- 
lation between the temperature and density is established 
(e.g. Theuns et al. 1998). The enhancement of QSO activi- 
ties peaking at 2 ~ 2 is expected to disturb this quiet state 
by the production of hard photons energetic enough to dou- 
bly ionize helium (e.g. Zheng et al. 2004). Modelling double 
helium reionization is the focus of the current paper. One of 
our basic assumptions is that QSOs are the main sources for 
He II reionization at low redshifts. This assumption is sus- 
tained by the evidence from observations pointing to He II 
reionization at redshifts concurrent with the escalation of 
QSO activities. Although star formation activity also peaks 
at similar redshifts, radiation from galaxies is too soft to 
contribute significantly to He II reionization (Benson et al. 
2002, Wyithe & Loeb 2003a), unless Population HI stars are 
considered (e.g. Venkatesan, Tumlinson & Shull 2003). 

In recent years commendable efforts have boon made to 
observe an He II Gunn-Peterson trough, at the Lya wave- 
length of 304A, in quasar absorption spectra. Estimation of 
the He II optical depth from the spectrum of the quasars 
Q0302-003 (2 = 3.286, Dobrzycki et al. 1991; Heap et al. 
2000) and of HE 2347-4342 (z = 2.885, Smcttc et al. 2002; 
Zheng et al. 2004) indicate a sharp decline in the optical 
depth at redshift z ~ 2.8. At 2 > 2.8 the optical depth r > 4 
and at 2 < 2.7 the optical depth r < 1. This decline implies 
that He II reionization occurred at 2 ~ 3. The optical depth 
of the quasar HS 1700-1-64 (2 = 2.743, Davidsen et al. 1996) 
T = 1.00 ± 0.07 at redshift 2 = 2.4, agrees with the general 
picture described above. It is important to stress that there 
may be other explanations for the large fluctuations in the 
He II optical depth. For example, in Smette et al. (2002), 
a model is described where optical depth fluctuations along 
the line of sight are a local effect due to neaxby "soft" and 
"hard" sources, and are unrelated to the global reionization 
of the IGM. 

We aim at a detailed modelling of the state of the IGM 
during He II reionization. Two approaches can be adopted 
in achieving this goal. The first is to employ hydrodynam- 
ical simulations that include radiative transfer and a good 
recipe for modelling the distribution of QSOs. An effort to- 
wards this goal has been made by Sokasian et al. (2001, 2002) 
who implemented their radiative transfer numerical code in 
a cosmological simulation of a box of 67h-iMpc (run by V. 
Springel). The simulations offer a valuable insight into the 



evolution of ionized regions and their topology. The effect 
of He II reionization on the thermal properties of the IGM 
have not been modelled in these simulations. This is be- 
cause radiative transfer is not merged self-consistently with 
the gas dynamics in the simulations. This is a technical lim- 
itation which is expected to be overcome in the near future. 
However, the simulations are limited by the available CPU 
power, preventing a thorough exploration of the key physi- 
cal parameters. The alternative approach is to develop fast 
semi-analytic models that capture the essential ingredients 
of the reionization process. This approach has been adopted 
to study a variety of other problems in structure formation, 
(e.g. Bi & Davidsen 1997; Kauffmann et al. 1997; Somcrville 
& Primack 1999; Kauffman et al. 1999; Cole et al. 2000; 
Benson et al. 2001). The basis of semi-analytic models is 
the synthesis of the knowledge acquired from diverse simu- 
lations and analytic reasoning. This synthesis allows us to 
study a given physical problem with greater detail than in- 
dividual simulations offer directly. In this paper we propose 
a Monte-Carlo model for studying He II reionization assum- 
ing QSOs are the sole sources of the ionizing radiation. For 
an assumed cosmology and a form for the QSO luminosity 
function, the model provides detailed information about the 
physical state of the IGM. The model is easy to run and al- 
lows a thorough exploration of the key physical process. We 
confront the model with available observational data, and 
provide predictions for the equation of state of the IGM, 
i.e., the local temperature-density relation (hereafter T — p 
relation). We compare the He II mean optical depth, t{z), 
to the optical depths calculated firom the He II Lya for- 
est of the three quasars mentioned above (Q0302-003, HE 
2347-4342 and HS 1700-1-64) and use this to constrain our 
model. We also compare the mean temperature, T{z), from 
our model to the temperatures at the mean density calcu- 
lated from nine quasars' Lya spectra (Schaye et al. 2000) to 
place further constraints. 

The paper is organized as follows. A short theoretical 

background is presented in § 2. In § 3 we describe our model 
which we use to investigate the process of patchy He II reion- 
ization and its effect on the thermal evolution of the IGM. 
In § 4 we try to fit observational measurements of the mean 
temperature and the mean optical depth of the IGM with 
our model for the power-law ACDM cosmological model im- 
phed by the WMAP data (Spergel et al. 2003). We also 
compute the local temperature-density relation for differ- 
ent redshifts. In § 4.2 we repeat the same computation and 
comparison for a ACDM model with running spectral index 
(RSI) (Spergel et al. 2003) and for an open CDM (OCDM) 
cosmological model. We conclude with a summary and a 
discussion of the results in § 5. 



2 AN OVERVIEW OF HE II REIONIZATION 

We assume that double helium reionization is mainly caused 
by ionizing photons emitted by QSOs, and that the con- 
tribution of galaxies is negligible. To justify this assump- 
tion we have computed the emission rates of He II ionizing 
photons from QSOs and galaxies and plotted them in fig- 
ure 1. The calculation of the emission rate of QSOs is based 
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on that of Madau, Haardt & Rees (1999; hereafter MHR) 
who performed a similar calculation for hydrogen (see Ap- 
pendix A). The emission rate of galaxies is estimated us- 
ing the GALFORM semi-analytic model of galaxy formation. 
Specifically, wc adopt the model and parameters of Benson 
et al. (2002), except for changes in the cosmological param- 
eters required to match the cosmological model considered 
here. The reader is referred to that paper, and references 
therein, for a full description of the model. We choose to ig- 
nore the cff'ccts of H I rcionization on later star formation* 
so as to obtain an upper limit on the number of ionizing 
photons produced at a given redshift. Note that this model 
tends to produce more photons than the more recent imple- 
mentation of the Durham model presented in Baugh ct al. 
(2004; see their Figure 1). As such, we are overestimating 
the number of ionizing photons produced. We assume that 
10% of the produced ionizing photons escape the galaxies 
into the IGM (observations of low and high redshift galax- 
ies imply escape fractions of < 10%; Leitherer et al. 1995; 
Tumhnson et al. 1999; Steidel, Pettini & Adelberger 2001). 
According to figure 1 the emission from QSOs is dominant 
at redshifts z > 1.4 even if only 20% of the QSOs emission 
calculated value reaches the IGM. 

At some high redshift QSOs begin ionizing He II. The 
fraction of He III prior to that redshift is assumed to be 
zero. The mean free path, /ph, of an He II ionizing photon is 



0.8n; 



jMpc where nue ii,-6 is the physical He II 



number density in 10~ cm~ . Initially the mean free path, 
/ph, is much smaller than the mean distance, dqso between 
QSOs and so only isolated patches (bubbles) of space are 
ionized. As the QSO activity increases more photons become 
available and the volume filling factor of ionized bubbles 
approaches unity. At this stage, the mean free path in the 
moderate density IGM, away from discrete absorbers, can 
exceed the Hubble radius, the radiation field becomes nearly 
uniform, and an equilibrium between recombinations and 
photo-ionization is established. This situation is maintained 
as long as enough photons are produced. 

When a given point in the IGM is ionized, its temper- 
ature increases by an amount proportional to the fraction 
of He II. This is a major source of heating in the IGM as 
each helium atom can be ionized and subsequently recom- 
binc several times (Miralda-Esc:udc & Rocs 1994). The dom- 
inant cooling mechanism is adiabatic cooling as a result of 
the cosmic expansion (e.g. Miralda-Escude & Rees 1994). 
During the patchy period of ionization large temperature 
differences between ionized and neutral regions can develop, 
spoiling any tight relation between the temperature and the 
density (T-p) in the IGM (e.g. Theuns et al. 1998). If this 
period is short lived adiabatic cooling and photo-heating in 
the later evolution of the IGM work to re-establish the tight 
T-p relation. Much of the thermal history of the IGM at 
2 < z < 5 is governed by the way He II reionization pro- 
ceeds. 
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* As shown by Benson et al. (2002), H i reionization leads to a 
suppression of star formation at later times. 



z - redshift 

Figure 1. The QSO and galaxy emission rate of He II ionizing 
photons per unit comoving volume. For QSOs we plot both the 
full value and 20% of the emission rate and for three different 
cosmologies. The solid line is for ACDM with Qq = 0.27, Ha = 
0.73 k. h = 0.72 (Spcrgcl ct al. 2003). The dashed line is for 
the running spectral index model of Spergel et al. (2003) with 
^0 = 0.268, = 0.732 h h = 0.71. The dotted line is for 
OCDM with Qo = 0.3, Ha = & h = 0.7. For galaxies we plot 
the emission rate with ionizing photons escape fra<;tion of 10% 
and running spectral index cosmological model. 



3 THE MONTE-CARLO MODEL 

Given N random points in the IGM with specified gas den- 
sity, temperature, and ionization state at some high redshift, 
«i, the model is designed to yield the corresponding quanti- 
ties at any lower redshift. The model assumes that a uniform 
ionizing background for hydrogen and singly ionized helium 
has already been established prior to Zi. The model aims at 
following the physical state of the IGM as a result of He II 
reionization by the QSOs population. At the early stages 
of He II reionization the mean free path for the absorption 
of an He II ionizing photon, Zph, is short and the reioniza- 
tion proceeds in a patchy way as a result of the discrete 
distribution of QSOs. At these stages a QSO generates an 
almost fully ionized bubble as the bounding thin (thickness 
Zph) ionization front propagates outward. The process con- 
tinues until the QSO fades away. The mean lifetime of a 
QSO is ~ lO'^ years (e.g. Hosokawa 2002; Schirber et al. 
2004; Porciani et al. 2004) which is much shorter than the 
recombination timescalc (order of 1 Gyr at « = 3) inside 
the bubbles. Therefore, we assume that the QSO radiation 
is emitted in bursts. At later stages when the filling factor of 
He III regions becomes significant, the mean free path, Zph, 
becomes so large that a non-negligible fraction of the ab- 
sorbed photons at a given time have actually been produced 
at significantly larger redshifts. We start the ionization of 
He II by QSOs at redshift Zi w 6, with hydrogen and he- 
lium already singly ionized. We neglect the He II ionization 
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Figure 2. Flowchart of the Monte-Ccirlo method. 

by galaxies. Figure 2 presents a schematic description of the 
Monte-Carlo model. The details of the model are as follows: 

(i) Initialization: we work with N points represent- 
ing random uncorrclatod points in the IGM (wo use A'^ = 
2 X 10^). At each point we determine the He III abundance 
fraction, X, the temperature, T, and the density perturba- 
tion, S. At the initial time, we set X = for all the points 
and assign temperatures according to T = To(l + S)'^^'^. At 
each point, the initial density contrast iSJ" at some very high 
redshift, z'" = 20, is drawn randomly from a normal (Gaus- 
sian) probabiHty distribution function with zero mean and 
with rms value ao. 

(ii) The local density and its evolution: the ionized 
fraction and the adiabatic cooling/heating depend on the 
local gas density and its evolution in time. Given the ini- 
tial mass power spectrum of the dark matter we would like 
to have a recipe for deriving the gas density contrast as a 
function of time. Hydrodynamical simulations strongly sup- 
port the picture in which the gas distribution in the IGM 
traces the dark matter density smoothed over the comov- 
ing Jeans length scale, xj, at which pressure forces roughly 
balance gravity (Hcrnquist ct al. 1996; Thcuns ct al. 1998). 
We therefore estimate the gas density power spectrum by 
smoothing that of the dark matter on the Jeans length scale 
(see Appendix C for details). We approximate the density 



contrast 5i{z) at a point i as follows. The first step is to 
obtain the gas density variance <|5?(«)^ = cr^{z) as a func- 
tion of redshift. Given the linear power spectrum we use 
the recipe described in Peacock (1999) to derive the corre- 
sponding nonlinear powcr-spcctrum. This power-spectrum 
is then multiplied by an appropriate Jeans smoothing win- 
dow and integrated to yield a^. Let be a sufficiently high 
redshift such that cr^(z'") <C 1 and such that the QSO emis- 
sion is negligible at z'". Draw a random variable 5p from 
a Gaussian distribution of zero mean and variance a (2'"). 
We write the density contrast, Si{z), at any redshift z < 
as 



il + 6i{z)) = e' 
where the nonlinear growth factor, D{z), is 



D{z) 



ln(l-hcT2(z)) 

ln(l-|-(T2(^in))- 



(1) 



(2) 



The form (1) is based on the log-normal model (e.g. Kof- 
man ct al. 1994) for the evolution of the density contrast, 
except that in the standard log-normal model the factor D 
is replaced with the linear growth factor of density perturba- 
tions. The form (2) for D ensures that the variance of 6i{z) 
is o-^{z) which can be estimated reliably using the recipe 
outlined below. On can show that (^5i{z)^ — <j^{z), as fol- 
lows. If (5p is Gaussian, it is easy to see that the mean of 
5i(z) is zero and the variance is exp \ji){zY (((^P)^)] — 1. 
Hence, using {{Sff) = cr^(«'") < 1 and (2), we find that 
(z)^ = cr^(z). Therefore, the form (1) guarantees a vari- 
ance of cr^(z) for Si{z), provided that 
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(iii) The ionization probability: at each time step the 
ionization state of helium is updated as follows. We assign 
a probability, fi, for photo- ionizing He II at the point i, and 
draw a random number, pi, between to 1 from a uniform 
distribution. If pi < fi then the point i is ionized. We adopt 
the following form for the photo-ionization probability. 



fi 



N Pbias(5i) 



(4) 



where riHc is the mean comoving number density of all he- 
lium species, Xj is the He III abundance fraction in point 
j, and Afabs(z) is the total number of ionizing photons per 
unit comoving volume that are actually absorbed in the cur- 
rent redshift interval [z{t), z{t + At)] (see below). At the be- 
ginning of each time step the number of available ionizing 
photons is computed as A^trs -I- NqAI where A'trs is the num- 
ber photons that have been transmitted from previous time 
steps and NqAt is the number of photons that are produced 
by QSO in the current time steps. The calculation of Nq is 
detailed in Appendix A. The number of available ionizing 
photons is tabulated in frequency bins and A^abs is computed 
in a self-consistent way as described in Appendix B. 

The factor Pbias takes account of the bias between the gas 
density distribution and the distribution of the He II ionizing 
photons. At the initial stages of He II reionization, the ioniz- 
ing photons are preferentially absorbed near the sources. As 
He II reionization proceeds, the mean free path, Zph, becomes 
large and the ionizing radiation eventually forms a uniform 
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background. In this case the distribution of the ionizing pho- 
tons is independent of 5. Here we adopt the following form 

for Pbias, 



M) « (l + <5.)''"''^ 



(5) 



where b is the bias parameter between the distribution of the 
QSOs and the gas density field, and s = s(Zph) is in general 
a function of the mean free path, Zph, of an He II ionizing 
photon. We work with three choices for s(Zph), 



Sconst — 1 



Slin 



and 



Xj/lph, if Iph > xj 
1, otherwise 



e(i-'phA.i)^ ifZph>a;j 
1, otherwise 



(6a) 



(6b) 



(6c) 



(iv) He II photo-ionization heating: once a point is 
flagged for ionization we update its ionized fraction and tem- 
perature according to, 

1, (7) 



Xi - 
and 



Ti+ATil-Xf""), 



(8) 



where AT is the increment in the temperature due to photo- 
ionization heating, and Xf'^"'^ is the ionization fraction be- 
fore the current ionization. We also make an attempt at 
modelling the effect of a finite QSO lifetime, tq, on the reion- 
ization process. The model does not contain explicit infor- 
mation on the spatial distribution of the sources. The way we 
alleviate this problem is by assuming that each freshly ion- 
ized point remains ionized for a period of time equal to Iq if 
it has not been selected for ionization during this time. Dur- 
ing this period of time, for each time step, some recombina- 
tion and photo-ionization heating occur, and can be treated 
as a local equilibrium. 

The process of drawing random numbers for points con- 
tinues until all A'abs photons are absorbed. 

(v) He II recombination: recombination reduces the 
He III abundance fraction by 

dX{l + 6) 



dt 



(9) 



where a^j^ jj is the He II case B recombination coefficient 
(e.g. Vomer & Ferland 1996) and rie in the mean comoving 
number density of electrons. 

At every time step we approximate the reduction in the 
He III abundance fraction at every point i due to recombi- 
nation 



Xi 



Xicxp -a'i^l „(T>)(1 + z)''(l + Si)n,At 



(10) 



(vi) Adiabatic cooling/heating: the main cooling pro- 
cess in the moderate density IGM at redshifts 2; < 6 is adi- 
abatic cooling. The adiabatic cooling/heating is calculated 
from the thermodynamic dependence between the tempera- 
ture and density of a monatomic ideal gas 
2/3 

(11) 



(Po) 



where To and (po) are the temperature and the mean density 
at some initial redshift, zo- The mean density, (p) ~ (l + z)^ 
so that 



(12) 



Therefore, from equation (12), the adiabatic cool- 
ing/heating at every time step is 



Ti{z) = Ti{z + Az)x 

( 1 + ^ Y( 1 + M^) ns) 

\l + z + AzJ \l + 6iiz + Az) J ' ^ ' 

(vii) Photo-ionization heating by H I &: He I: our 

basic assumption is that at redshifts z <Q hydrogen is fully 
ionized and helium is singly ionized. We also assume a local 
thermal equilibrium between the background radiation and 
recombination at every time step. For simplicity we treat 
He I as extra H I particles, and use the Miralda-Escude & 
Rees (1994) analysis for the IGM heating by the H I photo- 
ionization background radiation, to calculate the tempera- 
ture increment at every point i 

2 



Ti ^ Ti -h ^Ea^l\{Ti){l -h zf{l + Si)n,At, 



3k 



(14) 



where E is the average energy of the absorbed photons minus 
the average energy lost per recombination, and a^'j is the 
H I case B recombination coefficient (e.g. Verner & Ferland 
1996). 

(viii) Compton cooling: Compton cooling resulting 
from free electrons scattering off the CMB photons is non- 
negligible at redshifts 3 < « < 6. This cooling timescale is 
given by (Peebles 1968) 



ompton ■ 



1161.3(1 -FX"^) 



rGyr, 



(15) 



{l + zy[l-T^^^il + z)/T,] 

where Xe is the fraction of free electrons, jjp'^^ is the tem- 
perature of the CMB at the present day, and Tc is the tem- 
perature of the free electrons. At every time step the tem- 
perature decrease at every point i due to Compton cooling 
is 

T,At 



Compto 



(16) 



(ix) Other cooling processes: we have assessed the 
importance of the following cooling processes: coUisional 

ionization, recombination, dielectronic recombination, col- 
lisional excitation, and bremsstrahlung (e.g. Theuns et al. 
1998). We have found that all of these processes are negli- 
gible for S < 100 in the temperature range of interest to us 
T < lO^K. One should also keep in mind that we assume 
an instantaneous photo-ionization and photo-heating (equa- 
tion 7). Rapid energy loss in line excitations and coUisional 
ionization immediately follows and brings the temperature 
to T < lO^K. This rapid cooling is implicitly included in the 
temperature increment AT in equation (8). The parameter 
AT is treated as a free parameter in our model. 



3.1 The model input 

A basic input of the model are the cosmological parameters. 
We will present results for three variants of the Cold Dark 
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Table 1. Parameters of the three cosmological models considered 
here. The pajrameters of the power-law ACDM and the running 
spectral index ACDM models are taken from Spergel et al. (2003). 
The open CDM model has zero dark energy component. The pa- 
rameters Qq, f^A and Sib ^'I'c, respectively, the mean density (in 
units of the critical value) of the of the dark matter, the dark 
energy (cosmological constant) and the baryon component. The 
Hubble constant, h, in units of 100 km/s/Mpc, n is the power 
spectrum power-law, and erg is the rms value of density fluc- 
tuations in spheres of radius 8 h~^Mpc. The spectral index of 
the RSI model is a function of wavenumber, k, and is given by 
n{k) = n(fco) + [dn/dlnfe]ln(fe/A;o) with ko = 0.05 Mpc"'^ and 
dn/dlnfe = -0.031. 





no 


Qa 


fife 


h 


n 


0-8 


ACDM 

RSI 

OCDM 


0.270 
0.268 
0.300 


0.730 
0.732 
0.000 


0.0463 
0.0444 
0.0463 


0.72 
0.71 
0.70 


0.99 
Var. 
1.00 


0.90 
0.84 
0.90 



Table 2. The input parameters of the Monte-Carlo model: fj^^ 
is the QSO emission rate factor, tq is the mean QSO lifetime, 
To (2=6) is the initial mean temperature at redshift 2 = 6, AT is 
the temperature increment induced by photo-ionization of He II, 

E is the cross section weighted excess photon energy for hydro- 
gen ionization, xj{z=3) is the comoving Jeans length at redshift 
2 = 3, 6 is the QSOs distribution bias parameter, and s is the 
scaling function between the QSOs and the ionizing emission dis- 
tributions. 



Parameter 


Values 


Units 




0.2, 0.3, 0.4, 0.5, 1.0 




*Q 


0, 5, 10, 20, 40, 60 


Myr 


To (2=6) 


1.0, 2.0, 3.0, 4.0, 5.0, 6.0 


lO-^K 


AT 


1.0, 1.5, 1.84 


lO-'K 


E 


1.27, 2.54, 3.81 


eV 


a; J (2=3) 


0.05, 0.1, 0.2, 0.3, 0.4 


h-iMpc 


b 


1.0, 2.0, 3.0, 5.0, 8.0 




s 


const, linear, exp 





Matter (CDM) cosmological model. The first is the ACDM 
cosmological model of a flat universe with a cosmological 
constant and a power-law power spectrum of density fluc- 
tuations. The second is also a ACDM model with similar 
cosmological parameters but with a running spectral index 
(RSI) for the power spectrum. The third is an open CDM 
(OCDM) model, similar to our ACDM model but with no 
cosmological constant term. The main parameters deflning 
these cosmologies arc listed in Table 1. The parameters of 
the ACDM and RSI cosmologies are taken from the first 
year WMAP data analysis (Spergel et al. 2003). The OCDM 
model is not currently viable in view of the WMAP measure- 
ments and the observations of Type la SNe (e.g. Knop et al. 
2003; Ricss et al. 2004). Nevertheless, wo run our Montc- 
Carlo model with this cosmological model for the sake of 
comparison. 

The Monte-Carlo model has also to be fed with the 
input parameters pertaining to the physical processes dis- 



cussed previously. These parameters are summarized in Ta- 
ble 2. The first parameter in this table, fj^,^ , is the QSO 
emission factor. In Appendix A we describe our calculation 
of the QSO emission rate of ionizing photons per unit comov- 
ing volume, Nq{z). The calculation is based on the MHR lu- 
minosity function. Although we work with the basic shape of 
the inferred Nq (z) we introduce the factor /jy^ to allow for 
an overall reduction in the intensity of emitted radiation. 
The motivation behind this factor is that several observa- 
tions and simulations suggest a softer UV radiation field 
than the standard UV background due to QSOs emission 
(MHR). Smette et al. (2002) has shown regions of the HE 
2347-4342 spectrum in which the He II optical depth is very 
high and the corresponding H I optical depth is low required 
photo-ionization by a very soft UV radiation field. Indepen- 
dent evidence is provided by the column density ratios of 
highly ionized metals (A''si iv/Nc iv versus Nsi u/Nc iv) at 
high rcdshifts (Giroux & ShuU 1997; Savaglio et al. 1997; 
Songaila 1998; Kepner et al. 1999). Numerical simulations 
of the Lya forest also suggest a softer UV radiation field 
(Theuns et al. 1998; Efstathiou et al. 2000). 

Next is tq, the mean QSO lifetime. Since the points 
in the model are not spatially correlated we can take into 
account the QSOs lifetime only in a very limited way. We 
Eissume that a point that was ionized will stay ionized at 
least for the QSO lifetime. Therefore, = merely implies 
that the QSO lifetime is shorter then the model time step. 

The mean temperature at redshift z = 6, To(z=6), is 
fixed prior to He II reionization. This temperature is set by 

the interplay between cooling (mainly adiabatic and Comp- 
ton) and the photo-heating by the hydrogen ionizing back- 
ground. 

The parameter AT is the increment in the temperature 
due to He II reionization (cf. equation 8). In the limit of a 
very high ionization fraction an absorbed photon produce 
only one free electron which loses its energy by Coulomb 
interactions in the ionized IGM. The temperature increment 
in this case is given by (Miralda-Escude & Rees 1994) 

AT = (<e> - eo) ^ 1.84 x lO^K, (17) 

where X = nHe/{nn + nne + ne) w 0.07 is the fraction of 
particles which participate in the photo-ionization process 
from the total number of particles in the gas, <e> ~ 1.96 x 
10"^'^ erg is the mean energy per photon of photons with 
energies greater than 4 Ry which are emitted by QSOs, and 
£0=4 Ry ~ 8.72 x 10~^^ erg is the energy spent per He II 
ionization. 

However, when the ionization fraction is not very high, 
the electrons produced by photo-ionization may have a large 
probability of interacting with He 1 and lose energy in line 
excitations and collisional ionization. Furthermore, when 
secondary electrons are produced, the initial energy must 
be shared among several electrons, and so the final energy 
per electron is reduced. Therefore, we also explore lower val- 
ues of AT. 

The cross section weighted photon energy, E, deter- 
mines the photo-heating due to the H I and He I ionizing 
photon backgrounds, assuming local equilibrium between 
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the background radiation and the recombination processes 
(equation 14). Following Miralda-Escude & Rees (1994) we 
neglect the temperature dependence of E (which is valid for 
low temperatures, since the energy lost in recombination is 
then negligible). In the case of pure hydrogen and for QSOs 
SED L{iy) oc Miralda-Escude & Rees (1994) found 

E = 0.28 Ry = 3.81 eV. Since E should also include the He I 
ionizing photon background and since we use L{i') oc i/"^'* 
(see Appendix A2), we explore lower values of S as well. 

The next parameter, xj{z=3), is the Jeans length (co- 
moving) at redshift 2 = 3^. The Jeans length, xj, is de- 
fined as the scale over which the dark matter density should 
be smoothed to yield an estimate of the gas density (see 
Appendix C). Large values of xj yield smaller fluctuations, 
which in turn decreases the recombination rate, and also af- 
fects the adiabatic cooling. We assume a simple evolution of 
xj which depends only on redshift 

xj{z) = xj{z=3) j . (18) 

The QSOs bias parameter, b, determines the QSOs dis- 
tribution from the gas density distribution, and the scaling 
function, s, scales the ionizing emission distribution from 
the QSOs distribution. We used three different scaling func- 
tions constant, linear and exponential (see equations 6a, 6b 
& 6c). 



4 RESULTS 

We run our Monte-Carlo model for a wide range of the input 
parameters listed in Tabic 2. For every set of parameters, 
the Monte-Carlo model provides: 



Q0302-003 (Dobrzycki et al. 1991; Heap et al. 2000), and HE 
2347-4342 (Smette et al. 2002; Zheng et al. 2004). 

Unfortunately the available observational data are very 
limited. The uncertainties in the mean temperature data are 
large, while the uncertainties in the mean optical depth data 
are fairly small. There is a huge scatter in the optical depth 
measurements around redshift z = 2.8 (i.e. a scatter larger 
than expected on the basis of the reported error bars) and 
it is therefore not meaningful to use these measurements to 
compute a mean optical depth. Nevertheless the data still 
provide useful constraints on the model input parameters. 

As a measure of the "goodness of fit" of a given choice of 
parameters we use a statistic. For the temperature data 
we define Xt as 

2 X (^data ^modcl) /in\ 

XT = 4 ' (19) 

where the Tdata are the measured temperature in the data 
and Tmodei are the model predicted temperatures at the 
same redshifts as the data. The estimate of the errors in 
the data, ctt, (Schaye et al. 2000) should serve as a general 
guide only. We also define a x^ statistic for the optical depth 
measurement. We interpret these data as measurements of 
the local optical depth rather than the mean optical depth. 
Therefore, we write x? as 

2 X ^ ('^data Tmodel) 

^- = 1^ TP — rp — ) ' (20) 

where the Tdata are the measured optical depth in the data 

and Tjnodci are the model predicted mean optical depths at 
the same redshifts as the data. The observational errors in 
the data are cr r^ata the scatter about the mean optical 
depth as estimated in the model is crr^^odei — ((''' ~ ('''))^)^^^ 
(see figure 4) . 



(a) The He III actual volume filling factor and the volume 
filling factor of regions with He III fraction above 0.9 as a 
function of redshift. 

(b) The mean clumping factor for all regions and for re- 
gions with He III fraction above 0.9 as a function of redshift. 

(c) The mean temperature as a function of redshift for 
the mean density regions (regions with \S\ < 0.2) and for all 
regions. 

(d) The He II mean optical depth as a function of redshift. 

(e) The IGM equation of state, i.e., the T — p relation. 

We compare our results with available observational 
data. These data include (i) the mean temperature in re- 
gions of densities around the mean {\5\ < 0.2) estimated 
from nine quasars Lya spectra (Schaye et al. 2000), and [ii) 
the He II mean optical depth, tho ii(z) measured from the 
spectrum of the quasars HS 1700-1-64 (Davidsen et al. 1996), 

t Although the Jeans length is in principle specified by the tem- 
perature of the IGM we choose to treat it here as a free parameter 
of the model. The shape of the Jeans filtering window is quite 
uncertain — even in tlic huear regime, it depends non-trivially on 
the reionization history (e.g. Hui & Gnedin 1997, Nusser 2000). 
The situation becomes more complicated when nonlinear effects 
become important. 



4.1 Results with a power-law ACDM cosmology 

We first present results for the ACDM cosmological model 
(see Table 1). We discuss other cosmologies in section § 4.2. 

In general, the best fits of the model to the observations, 
with Xt < 1 and x? < li are obtained for /j^^ « 0.3, Jeans 
scale of xi{z = S) fs 0.2 - 0.4 h'^Mpc, and QSOs lifetime 
^ 10 Myr. We found no constrains on the bias parameter 
b, and the scaling function s. 

The initial mean temperature is To[z = Q) < 40000K. 
The temperature increment is AT Pi 10000- 15000K, which, 
as expected, is smaller than the estimate of Miralda-Escude 
& Rees (1994) who find AT ~ 18400K for a fully ionized 
regions. This is due to energy losses by line excitation and 
coUisional ionization in partially ionized regions. We found 
a correlation between the temperature increment, AT, and 
the hydrogen photo-ionization heating parameter, E. For 
lower temperature increment, AT « lOOOOK, the heating 
parameter is 1.27 ^ E < 3.81 eV, but for higher temperature 
increment, AT « 15000K, the heating parameter is only 
E fs 1.27 eV. 

Figure 3 summarizes some of our main results for the 
ACDM cosmology with parameters (listed in the top-left 
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Figure 3. Results for the power-law ACDM cosmological model for a choice of set of parameters as indicated in the figure. This set of 
parameters gives = 0.76 and Xt — 0.75. Panel (a): The He III fraction (see text) as a function of redshift. The solid line is the mean 
fraction of He HI. The dashed line is the volume filling factor of regions with He III fraction above 0.9. Panel (b): The mean clumping 
factor, (1 + a^), as a function of redshift. The solid line is the mean clumping factor. The dashed line is the clumping factor for regions 
with He HI fraction above 0.9. Panel (c): The He II mean optical depth as a function of redshift. The solid line is the mean optical depth 
calculated from the model. The symbols are measurements of the He II optical depth from the spectrum of the He II Gunn-Peterson 
effect towards the quasars HS 1700+64 (triangles, Davidsen et al. 1996), Q0302-003 (squares, Dobrzycki et al. 1991; Heap et al. 2000), 
and HE 2347-4342 (circles, Smettc et al. 2002; Zheng et al. 2004). The noise in the line is due to the limited number of points used in 
the Monte-Carlo model. Panel (d): The mean temperature as a function of redshift. The solid line is the mean temperature of the mean 
density regions, calculated from regions with |5| < 0.2. The dashed line is the overall mean temperature, calculated from all regions. The 
circles with error bars are the observed mean temperature calculated from nine quasars Lya spectra (Schaye et al. 2000). 



panel) for which the Monte-Carlo model yields one of the 
best fits to the observational data, Xr = 0-76 and x? = 0.75. 

Wc quantify the development of He III regions by com- 
puting the mean ionized fraction, X (shown as the solid 
curve in panel (a) of figure 3), and the fraction of volume 
occupied by regions having X > 0.9 (dashed curve). The 
dashed and solid curves behave as expected. The redshift at 
which the dashed curve reaches unity in panel (a) is z = 2.9. 

Panel (b) of figure 3 shows the evolution of the clump- 
ing factor, /clump = (l + a^), for all regions (solid curve) and 



for regions with X > 0.9 (dashed curve) as a function of red- 
shift. The clumping factor for all regions increases in time as 
expected. The clumping factor indicated by the dashed curve 
is large at high redshifts since the reionization starts in dense 
regions, and decreases in time until it meets the solid curve 
around the redshift where full reionization is established. At 
later times all the regions are almost completely ionized, and 
therefore the dashed curve follows the solid curve. 

Next we examine the optical depth, thb n, for helium 
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Figure 4. The rms values of optical depth and the temperature fluctuations in the power-law ACDM cosmological model for the set of 
parameters indicated in the figure. Panel (a): The He II optical depth rms as a function of redshift. The noise in the curve is due to the 

limited number of points used in the Monte-Carlo model. Panel (b); The rms of temperature fluctuations as a function of redshift. The 
solid line is the temperature rms in regions with |i5| < 0.2. The dashed line is the overall temperature rms, calculated from all regions. 
The rms value of the fluctuations in a variable Z is Zrms = {iZ — (Z))^)^/^. 
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Figure 5. The mean temperature as a function of redshift in the power-law ACDM cosmological model for two sets of parameters as 
indicated in the panels. The rest of the model parameters are the same for both cases: /jy^ = 0.3, tq = Myr, a;j(2=3) = 0.3, 6 = 5 & 

s = linear. The set of parameters of panel (a) gives = 0.93 and x't — 0.77. The set of parameters of panel (b) gives = 1-0 and 
Xr = 0.94. The solid line is the mean temperature of the mean density regions, calculated from regions with |<5| < 0.2. The dashed line 
is the overall mean temperature, calculated from all regions. The circles with error bars are the observed mean temperature calculated 
from nine quasars Lya spectra (Schaye et al. 2000). 



Lya absorption. This is given by 



(21) 



where H{z) is the Hubble constant at redshift z, and ai — 
ne^ /nieC^jXifij ^ 1.1 x 10~^* cm^ is the integrated Lya 
scattering cross section of the He II transition at wavelength 
A( = 304 A and oscillator strength fij . At every time step. 



we calculate the optical depth at each point i, and the mean 
optical depth 



THe 11 = — In 



1 / He II\ 

-2^exp(-r, ) 



(22) 



The solid line in panel (c) of figure 3 is THe u computed from 
our model as a function of redshift according to (22). The 
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Figure 6. The tcrnpcraturc-dcnsity relation for the power-law ACDM eosmological model. The values of the key parameters are indicated 
in the figure. The six panels show the temperature as a function of the density contrast S at different redshifts. The dots are the local 
temperature values. T^g is the effective equation of state assuming T^ff oc (1 + 5)^^^ for all points ionized and non-ionized. The effective 
temperature is calculated from the least squares linear fit between log(T) and log(l + <5). 



noise in the line is due to the limited number of points used 
in the Monte-Carlo model. For comparison we also show the 
available observational data points measured from the spec- 
trum of the He II Gunn- Peterson effect towards the quasars 
HS 1700+64 (triangles, Davidsen ct al. 1996), Q0302-003 
(squares, Dobrzycki et al. 1991; ffeap et al. 2000), and HE 
2347-4342 (circles, Smette et al. 2002; Zheng et al. 2004). 
One can easily see a significant decrement in the optical 
depth from t ~ 10 at z > 3.5 to t < 0.5 at « < 2.5. 

In panel (a) of figure 4 we see that the rms of the fluc- 
tuations of the optical depth, T^ms = ((t— (t))'^)^^^ , behaves 
in a similar manner to the mean optical depth. At high red- 
shifts, z > 3.5, the high density regions get ionized while the 
low density regions stay neutral. This results in a large scat- 
ter of the optical depth (Trms ^ 20). The scatter decreases 

towards full reionization^ . By 2 ~ 3, most of the high den- 
sity regions are already ionized as well as a large fraction of 
the low density regions. At ^; < 2.5 we have Trms ^ 1- 

Panel (d) of figure 3 shows the evolution of the mean 

■t By full reionization we mean a state in which the volume filling 
factor of regions with X > 0.9 is unity. 



IGM temperature as a function of redshift. The model pre- 
dictions are shown as the solid and dashed curves. In each 
panel, the solid curve refers to the mean temperature at 
points with density contrast \5\ < 0.2, while the dashed 
line is the volume weighted mean temperature defined as 
the average temperature of all the points of the Monte- 
Carlo model. Overlaid as the points with error bars are the 
temperatures as inferred from the observations (Schaye et 
al. 2000). The data points correspond to temperatures in 
regions of mean density, and therefore they should be di- 
rectly compared with the solid curves. However, the differ- 
ence between the solid and dashed curves in each panel is 
within the error baxs. The gradual variation of the model 
temperatures with redshift is in reasonable agreement with 
the (lata. Comparing with panel (a), we see that the tem- 
perature curves peak at the redshift for which full reioniza- 
tion is reached. At higher redshifts each He III particle can 
recombine more than once, resulting in an efficient photo- 
heating by the He II ionizing photons. At lower redshifts, 
the recombination is greatly reduced and, subsequently, the 
photo-heating rate decreases and adiabatic cooling domi- 
nates. At these low redshifts, the decline of the solid lines 
does not seem to be fast enough to match the observational 
data. The dashed curves, however, fall nicely among the data 
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Figure 7. The temperature-density relation as derived from two sets of observed absorption-line systems (Bryan & Machacek 2000) 
and from the Monte-Carlo model with the same parameters set as in figure 6. The effective equation of state, Tgff oc (1 + 5)^"^, is 
calculated from the least squares linear fit between log(T) and log(l + <5). Right panel: HS 1946+7658 at (z) = 2.7 (red circles), data 
from the Monte-Carlo model at the same redshift (black dots). Left panel: APM 08279+5255 at (z) = 3.4 (red circles), data from the 
Monte-Carlo model at the same redshift (black dots). 



points. This difference between the dashed and solid curves 
in each panel is due to the slower adiabatic cooling of re- 
gions with \5\ < 0.2, relative to the overall volume weighted 
cooling which is dominated by the expanding low density 
regions. 

Unsurprisingly we found a tight relation between the 
shape of the mean temperature curve and the combination 
of the heating sources AT and E. In the case we have investi- 
gated so far, with low temperature increment AT — lOOOOK 
and moderate hydrogen photo-ionization heating parameter 
E — 2.54 eV, the temperature rises in 5000 — 6000K during 
3 ^ z <j 4:, and never rises above 18200K. We also present 
two other typical cases in figure 5. In panel (a) of figure 5 we 
present a case of low AT = lOOOOK and high E = 3.81 eV. 
In this case the rise in the temperature during reionization 
is around 7500K. The high hydrogen photo-ionization heat- 
ing parameter does not allow the IGM to cool at higher 
redshifts, 4 > z > 6, and the temperature reaches a maxi- 
mum around 19000K. In panel (b) of figure 5 we present a 
case of high AT = 15000K and low E = 1.27 eV. In this 
case the rise in the temperature during reionization exceeds 
8000K and can reach llOOOK. The temperature can reach a 
maximum of 20000K. 

It is instructive to examine the temperature fiuctua- 
tions in space. This can be quantified in our Monte-Carlo 
model by computing the rms of the temperature fiuctua- 
tions, Trms = ((T - (r))^)^/^ In panel (b) of figure 4 we 
show the temperature rms as a function of redshift. The tem- 
perature rms of regions with mean density (the solid lines) 
is affected mainly by helium photo-heating. When the IGM 
reaches full ionization, the photo-heating is not effective any 
more and the temperature rms decreases to approximately 
the same value as before ionization. 



Finally, we consider the effect of He II ionization on 
the T ~ p relation. In the absence of He 11 reionization, we 
expect a tight T — p relation resulting from adiabatic cooling 
and photo- heating (e.g. Hui & Gnedin 1997). We assume in 
the Monte-Carlo model that at z = 6 the T — p relation is 
such that T oc (1 + 5)'^^^ . To emphasize the effect of patchy 
helium reionization we do not include any scatter in the 
T — p relation at z = 6. 

Figure 6 presents the results for the same set of pa- 
rameters as the above. The six panels in the figure, corre- 
spond to the T — p relation at different redshifts. At the 
bottom of each panel, an effective equation of state of the 
form Toff oc (1 + 5)^"^ is given. This effective equation of 
state is derived by the linear least square fit between log(r) 
and log(l + 5). 

In figure 6 we see that during the early stages of reion- 
ization, when mostly high density regions are being ion- 
ized, the slope of temperature-density relation is steepened 
and the scatter increases significantly. At z ~ 3 — 4 a non- 
negligible scatter of about 50% is introduced. Around z — 3, 
where the reionization of low density regions is much more 
significant, the slope of temperature-density relation fiat- 
tens. When the helium in the IGM is fully ionized the photo- 
ionization heating becomes inefficient and the main process 
affecting the gas temperature is adiabatic cooling. The cool- 
ing reduces the scatter, and by z = 1 a tight T — p relation 
is restored. 

In figure 7 we present the T — p relation as derived from 
two sets of observed quasar absorption-line systems from 
Bryan & Machacek (2000): APM 08279+5255 at (z) = 3.4 
(left panel) and HS 1946+7658 at (z) = 2.7 (right panel), 
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Figure 8. The slope, 7, of the eflfective equation of state, Tgn oc 
(1 + 6)1-'^, at the mean density. The open circles are the slopes 
calculated from the nine QSO Lya spectra in panel (d) of figure 3 
by Schaye et al. 2000. The filled circles are the slopes calculated 
from the Montc-Carlo model for regions with |5| < 0.2 at the 
same redshifts. la error bars in 7 are presented for both the ob- 
servational data and the model points. 



Table 3. Comparison between the slopes, 7, of the effective equa- 
tion of state at the mean density from the observations (Schaye 
et al. 2000) and the slopes from the Monte-Carlo model (from 
regions with |5| < 0.2). 2 is the redshift with the errors from the 

observational data. 7data is the slope from the observational data. 
7modei is the slope from the Monte-Carlo model, rj is the distance 
in a between the model slope and the observational data slope. 
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and compare with the T — p relation from the Monte-Carlo 
model at the same redshifts. 

Ax z = 2.7 the slope of the effective equation of state 
from the model, 7modei = 1.392 ± 0.002, fits the slope of 
the data from the HS 1946+7658 absorption-lines, 7data = 
1.31 ±0.15, within O.Gcr. At z = 3.4 the slope of the effective 
equation of state from the model, 7modoi = 1.748±0.008, fits 
the slope of the data from the APM 08279±5255 absorption- 
lines, 7data = 1.55 ±0.14, within 1.4(7. The coefficient of the 
effective equation of state from the model fits the data from 
the observations within 1.2ct sX z = 2.7 and within 0.3a" at 
z = 3.4. 

In figure 8 we present the slopes, 7, of the effective equa- 
tion of state at the mean density of the nine quasar Lya 
spectra studied by from Schaye et al. (2000; open circles) 
and the slopes calculated from the Monte-Carlo model from 
regions with \S\ < 0.2 at the same redshifts (filled circles). As 
expected, the reheating of the gas at the time of reionization, 
2: ~ 3 — 4, causes a decrease in 7. After the reionization the 
IGM gradually evolves again towards an ionization equilib- 
rium and 7 increases. To compare between the slopes from 
the observational data and the slopes from the model we 
calculated 77, the distance in a between the two 7, 

_ lldata - 

^ ~ I 1 2 ' ^ ' 

V "^data + '''model 

where 7data and iTdata arc the observational data slopes and 
their Icr errors, and 7modei and (Jmodei are the Monte-Carlo 
model slopes and their Icr errors (see Table 3). Prom figure 8 
one can see that in general 7 from the Monte-Carlo model 



have higher values then 7 from the observational data, al- 
though for most redshifts 77 < 3 (see Table 3). 



4.2 Results with the RSI and OCDM cosmologies 

We have also run our Montc-Carlo model for the RSI and 
the OCDM cosmological models (see Table 1) and have tried 
to fit the model parameters to the observations of the mean 
optical depth and mean temperature. We found no major 
differences among the results in the three cosmological mod- 
els we consider 

As in the ACDM cosmology, by fitting the observations 
with Xt ^ 1 and Xt ^ 1, we get /j^^ « 0.3 for both the 
RSI and the OCDM models. The QSOs lifetime in the RSI 
and the OCDM models is < 5 Myr. A Jeans scale length 
of a;j(2 = 3) ^ 0.1 - 0.4 h^^Mpc fits the RSI model, while 
2:j(z=3) ?a 0.2-0.4 h^^Mpc fits the ACDM, and only a;j(z= 
3) « 0.4 h-^Mpc fits the OCDM model. While in ACDM 
the bias parameter is not constrained, we find 6 « 2 — 8 and 
6 « 2 — 3 for the RSI and OCDM models, respectively. As in 
the ACDM cosmology, we found no constrains on the scaling 
function (equation 6). In the RSI and OCDM cosmological 
models the initial mean temperature is To(a=6) < 30000K, 
and temperature increment of AT « 10000 - 15000K. For 
AT ~ lOOOOK the hydrogen photo-ionization heating pa- 
rameter varies between 2.54 < S < 3.81 eV while in the 
ACDM model 1.27 <E< 3.81 eV. In all cosmological mod- 
els, the heating parameter is 1.27 eV for AT « 15000K . 
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Figure 9. The filling factor, temperature and optical depth evolution as a function of redshift in a ACDM cosmology for / 



(right panels) /j^ 



0.2 (left panels). For /^^ 



1.0 the set of parameters gives Xt ~ 2-3 and Xt — '^■^t ^'^d for 



: 0.2, Xt ■ 



1.0 
: 3.2 



and Xr 



: 0.1. 



5 DISCUSSION 

We have presented a model for following the thermal proper- 
ties of the IGM during patchy He II roionization. The model 
assumes that radiation from QSOs is the main cause of 
He II reionization and neglects the contribution from galax- 
ies. This is consistent with our finding that the model can 
account for He II reionization in a way consistent with the 
available observations of the mean temperature and He II 
optical depth. Further, according to the GALFORM semi- 
analytic model for galaxy formation the contribution of 
gala:x;ies to the He II ionizing radiation is negligible com- 
pared to QSOs even at redshifts as high as ~ 6. We have 



also neglected contributions from thermal emission in shock 
heated gas (Miniati et al. 2004) and also from redshifted 
soft X-ray produced at 2 ~ 20 (Ricotti & Ostriker 2004). 
These contributions could be important for He II, especially 
at 2 > 3, and will be studied in detail in a future paper. 

We have considered three different cosmological models: 
standard power-law ACDM, ACDM with a running spectral 
index (RSI), and an open CDM and found no significant dif- 
ferences between them. Table 4 lists the input parameters 
for which the Monte-Carlo model matches best the observar 
tional data. 

A key parameter is the emission factor fj^^ which mul- 
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Table 4. Summary of the best fit parameters derived from com- 
parison with the observed optical depth and temperature evolu- 
tion for Xt < 1 and Xt ^ 1- "N.C." implies no constraints on the 
parameter in the relevant range. 
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tiplies the overall emission rate derived from the MHR lu- 
minosity function. In all cosmologies, an emission factor of 
fj^^ = 1 yields results that are inconsistent with the rele- 
vant observational data (see panels on the right of figure 9). 
For /jy^ = 1, reionization occurs too early at z = 4.1, and 
shifts the typical behaviour of the mean temperature and 
optical depth to earlier times. Too low a value for /^^^ does 
not match the data either. For = 0.2 the reionization 
occurs too late (z « 2.5), and the mean optical depth does 
not reach high enough values around z — 3 (see panels on 
the left of figure 9). To match the data we need /^^ ^ 0.3. 
This result agrees with several observations (Giroux & ShuU 
1997; Savaglio et al. 1997; Songaila 1998; Kepner et al. 1999; 
Smette et al. 2002) and simulations (Theuns et al. 1998; Efs- 
tathiou et al. 2000) which suggest a softer UV radiation field 
than the standard UV background due to QSOs emission. 

The He III fraction rises modestly at the initial stages 
of reionization. The rise becomes more rapid as the He III 
fraction gets close to unity (full reionization), at redshift 
Zion. In the simulations of Sokasian et al. (2002) the full 
reionization occurs at a redshift of z « 3.8 (see their figure 
3, models 1-3), while we found that full reionization is never 
achieved before Zion = 3. The reason for this difference is 
that we tune our model to match the temperature evolution. 
In our model the temperature increases until full reionization 
is reached, and declines monotonically afterwards. 

The optical depth, for He II Lya absorption, is thb n ~ 
10 at a > 3.5 with a slow decline. Around Zion the decline is 
much more rapid, from thb ii ~ 10 at « ~ 3.5, to Tae ii ~ 0.5 
at 2 ~ 2.5. At z < 2.5 the decline is gradual again and 
reached thc ii ~ 0.3 at z — 1.0. We need more observational 
data to have better constrains on Monte-Carlo model param- 
eters. Sokasian et al. (2002) also compute the mean optical 
depth in their simulations, however, they show results up 
to 2 = 2.8. Some of their models also match the observa- 
tions up to that redshift, but we cannot make a more direct 
comparison of their models with ours at lower redshifts. 

An important application of the model is to study 
the temperature-density relation in the IGM. Patchy He II 
reionization is expected to produce a 50% scatter in this re- 



lation at 3 < 2 < 4 (sec also Hui & Haiman 2003). At lower 
redshifts, when the helium in the IGM becomes fully ion- 
ized, photo-heating is not effective anymore and the main 
process governing the gas temperature is adiabatic cooling. 
The scatter is then reduced and a tight T — p is rapidly es- 
tablished. The increase in the scatter at 3 < 2 < 4 should be 
taken into account in the analysis of the Lya forest data and 
can be detected in high resolution spectra. The temperature- 
density relation from our Monte-Carlo model fits the data 
from Bryan & Machacck (2000) around 2 ~ 3 reasonably 
well, but more observational data is needed, especially in 
low density regions, to be able to put better constrains on 
the model. 

The model predicts a gradual rise in the mean temper- 
ature between 2 ~ 4 and 2 ~ 3. Our results match the ob- 
served mean temperature as a function of redshift (Schaye et 

al. 2000). A possible discrepancy is that adiabatic cooling^ 
is not efficient enough to cause a rapid decline in the tem- 
perature at 2 < 2, as the observations indicate. However, 
the observational situation is still inconclusive. Although the 
temperature increase inferred from the observations is prob- 
ably robust (Theuns et al. 2002), the quantitative behaviour 
is less certain. 

There is a weak degeneracy between the temperature 
increment, AT, and thc energy, E, corresponding to photo- 
heating by the diffuse hydrogen ionizing background. An 
energy of i? ~ 1.27 eV works well with the full range 
10000 < Ar < 15000K. On the other hand the higher ener- 
gies E ~ 2.54-3.81 eV works only with AT close to lOOOOK. 
In any case, the inferred values for AT are consistent with 
the assumption that due to energy loses in partially ionized 
regions the temperature increment must be smaller then the 
estimates of MiraldarEscude & Rees (1994) for fully ionized 
regions. 

We have also found that only short QSO lifetimes fit the 
data: Iq < 10 Myr for the ACDM model and tq < 5 Myr for 
the RSI and OCDM models. The initial mean temperature 
has an upper limit of To (2 = 6) 40000K for the ACDM 
model and an upper Umit of To (2=6) w 30000K for the RSI 
and OCDM models. The Jeans scale length of a;,j(2=3) ~ 
0.1 - 0.4 h^^Mpc fits thc RSI model, while x,j(2=3) ^ 0.2 - 
0.4 h-^Mpc fits the ACDM, and only xj(2=3) « 0.4 h'^Mpc 
fits the OCDM model. While in ACDM the bias parameter 
is not constrained, in RSI & « 2 — 8 and in OCDM 6 « 2 — 3. 
In all cosmologies we found no constrains on the scaling 
function. 

In numerical hydrodynamical simulations of the Lya 
forest He II reionization is often modelled as a result of a uni- 
form ionizing background of photons. This recipe for He II 
reionization leads to a sudden jump in the mean IGM tem- 
perature at about 2 ~ 3. Our findings imply that this recipe 
is oversimplified and may lead to incorrect conclusions. An 
incorporation of models similar to ours in this type of sim- 
ulation is needed. 

The distribution of the column densities ratio r] = 



§ Other cooling mechanisms are inefficient at the relevant range 
of densities and redshifts (Theuns et al. 1998). 
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iV(He II)/A'^(H I) contains important information on the 
properties of the ionizing sources (e.g. FUSE observations, 
Zheng et al. 2004). Therefore, it will be very interesting to 
look at the distribution of r; as a diagnostic of the reioniza- 
tion epoch. To do so one must include H I rcionization, and 
therefore to take into account the ionizing radiation from 
galaxies. It is not trivial to add it to our Monte-Caxlo model 
since the model is based on the assumption that the ionizing 
sources are short lived, which is incorrect for galaxies. There- 
fore, adding the ionizing radiation from galaxies requires a 
different approach which will be presented in a future paper. 
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Table Al. Parameters of the double power-law luminosity func- 
tion (Pei 1995; MHR). 



r^iraiiK'lc]; X'ahic 

Pl 1.64 

132 3.52 

1.9 

C 2.58 

C 3.16 

M.(0) -22.35 

log[^,(h5o)/Gpc-3] ... 2.95 



APPENDIX A: THE QSO IONIZING PHOTON 
EMISSION RATE 

Our calculation of the QSO emission rate of Ho II ionizing 
photons per unit comoving volume is based on MHR for 
hydrogen. First we calculate the QSO emission rate in a flat 
universe without a cosmological constant (f2o = 1, f^A = 0). 
Later in Appendix A6 wc extend this calculation to fit all 
cosmologies with any desirable Oq and Qa. 



Al QSO luminosity function 

Following MHR we represent the quasax blue luminosity 
function (Ab = 4400 A) as a double power-law: 



<t>{L,z) 



(Al) 



where /3i and (52 are the power-law indices for faint and 
bright quasars, respectively. The position of the break L* {z) 
is 



L,{z) = L»{Q){1 + z) 
where 



(A2) 



L.(0) =L0XlO-°-*'^*(°)-^®l =4.67x10^^ erg Hz-^s'i, 

(A3) 

Mq = 5.48 and L© = 3.44 x 10^® erg Hz"^s"^ are the Solar 
magnitude and luminosity respectively, and the rest of the 
luminosity parameters can be found in Table Al. 

In Table Al is given for a Hubble constant h = 0.5. 
From this value we calculate 0* (ft) as a function of h: 



,{h) = 8.9125 X (^Aj X 10"'' Mpc" 



(A4) 



L(i') oc v , with different slopes in different wavelength 
ranges (MHR), 



L(i.) 



(2500 A < A < 4400 A) 
(1050 A < A < 2500 A) 
(A < 1050 A) 



(A5) 



The QSO ionizing flux, S, is the number of ionizing photons 
emitted by a QSO per second. 



hv' 



(A6) 



where h is the Planck constant. For convenience wc represent 
5* as a function of wavelength, A, and since the ionization 
wavelength of helium (A = 228 A) is less then 1050 A, 



/ 4400 \ - 



/2500\-°* / 



1050 \ ■ 



1.8/1 V250oj U05oj V A j 



228 A 

oA (A7) 



where Ls is the luminosity in the blue band (in units of 
erg Hz~^s~^), and neglecting ionization of heavier elements. 



A3 QSO ionizing flux function 

Equations (Al) and (A7) allow us to write the number of 
emitted ionizing photons per unit time per unit (comoving) 
volume as 



(l>*/S4z) 



where S* {z) for helium ionization is 

L«(z) /4400\-°-» /2500\-"« /1050\-i' 

~ Tsft" I2500; lio5o; \~) 



= 1.057x10^^(1-1-0)' 



phot s ^ . 



(AS) 

228 A 

A 

'(A9) 



A2 The QSO spectral energy distribution 



A4 K-Correction 



The luminosity spectrum of a "typical" QSO is assumed 
to have a power- law spectral energy distribution (SED), 



Taking into account the wavelength dependence of the 
power as in the QSO spectrum, one should replace the K- 
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correction (1 + 2)"° ^ in equation (A9), with 

(1 + zf'-"-^ (2500 A < i^SoA < 4400 A) 

1+ 



(1 + (1050 A < i^23A < 2500 A) , 

(AID) 



. KiK2il + zy-^-^ (i 



l + z 



< 



1050 A) 



where if, = (|i§)« -°-%nd;^. = 



0.8-1.. 



A5 The QSOs emission rate of ionizing photons 
per unit comoving volume 

The QSO emission rate of ionizing photons per unit comov- 
ing volume can be written as 



Nq{z)= I ^{S,z)SAS. 



(All) 



where Smin = 0.0185* (0) corresponds to Mb ^ —18, and 
therefore L-mir, ~ 0.018L,(0) (Cheng et al. 1985). Using 
equation (AS) we write Nq{z) as 



Jo. 018 ^ 

where x = S/St{z). Numerical integration gives 

xdx 



: 2.298. 



(A12) 



(A13) 



Equations (A12), (A4) and (A9) yield the QSO emission 
rate of helium-ionizing photons 

iV,(.) = 2.16xl0-X,(.)^^±^(A)%-Mpc-. 

e? +e'= V0.5/ ^^^^^ 



A6 The cosmological factor 

Untill now we assumed a flat universe cosmology with 
f2o = 1 and 57a = 0. To obtain Nq{z) in different cosmolo- 
gies we multiply the above result by a cosmological factor 
fco3mo{^o,^A) so that 



Nq{z) =^ fcosmo{iio,^A)NQ{z), 

where 



(A15) 



/cosmo(^0? ^a) — 



drc(l,0)/dg 
drc(no, QA)/dz 

-\/(l — a)Qo + {a^ — a)QA + a, 



(A16) 

and drc/dz is the rate of change of comoving distance with 
redshift. 



APPENDIX B: THE NUMBER OF ABSORBED 
PHOTONS PER REDSHIFT INTERVAL BIN 

Wc write the total number of ionizing photons that were 
absorbed at time interval [t,t + At] as 

A^abs(«) = [NQ{z)At + Nu.] [1 - e-^^^/'pi-f-)] , (Bl) 



where Nq{z) is the QSO ionizing photon emission rate per 
unit comoving volume, At is the time step, A'trs is the num- 
ber of photons received from previous time intervals, and 
the mean free path of the He II ionizing photons is 

/ph(a) = ["He Il(l + 5{z)){l + zfaion] , (B2) 



where nne 11 is the mean comoving number density of He II, 
5{z) is the local density perturbation, and CTion is the mean 
ion-photon cross section for He II. 

The ion-photon cross section at frequency u is 



,(i/) = 8.61 X 10"^® cm^ (—^) 



-4 4arctane/e 

1 - ' 



(B3) 



where e = [I'/i^e n — 1] and me n = 1.3 x^® Hz is the 
frequency at the He II ionization energy, E = 54.42 eV (Cen 
1992). 

Following Theuns et al. (1998) we approximate the Cen 
(1992) ion-photon cross section as having simple power- low 
dependence on frequency 



^{u) = 1.9 X 10"^® cm^ 



\ t^e II / 



(B4) 



Bl Fi"equency binning 

Because of the frequency dependency of the cross section 

we divide the frequency interval of ionizing photons into 10 
bins, from i^hc ii to infinity, and calculate the total number 
of absorbed photons in each bin. Wc define the luminos- 
ity as the photons energy rate per frequency (in units of 
erg Hz~^s~^). The luminosity of a "typical" QSO at wave- 
lengths A < 1050 A is assumed to behave as L(u) oc u~^'^ 
(equation A5). To have equal luminosity in each bin, we 
choose the frequency at the edge of every bin j to be 



y = 1^1050 



I'lOBO 



( 



V I'lOSO 



-) 

/ 



-1/1.; 



(B5 



where i^ioso = 2.85 x lO^^Hz is the frequency at wavelength 

A = 1050 A, and /b = 0.1 is 1 over the number of bins. For 
each bin we calculated the mean ion-photon cross section 



ion-'''^ VHe II / Vi'lOSO/ 



J 1/3 V I'lOSO-' 



O-o (fHe 11) — 7+T 



^dl^ 



(B6) 



(B7) 



and the QSO ionizing photon emission rate 
B2 The number of absorbed photons 



To get more accurate results for the radiative transfer we 
first calculate the number of redshifted photons that were 
transferred to a lower frequency bin. Next, since we choose 
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10 frequency bin we solve a set of 10 differential equations 
for the photon transmission, 



10 



V j=i ) (B8) 

where A/'o(«) = ^ N^{z)At+Ni^^{z) is the total number 
of photons that are available for ionization and N^^^ is the 
previous number of transmitted photons in bin j. The actual 
number of the absorbed ionizing photons is 



iVabsW=iVo(^)-^iV4(^). (B9) 



APPENDIX C: THE GAS DENSITY FIELD 

We assume that the gas density traces the dark matter den- 
sity smoothed over the Jeans length scale, xj, determined by 
the balance of pressure and gravitational forces. Under this 
assumption wo compute the variance of the gas density field, 
cr^, from the dark matter nonlinear (dimensionless) power 
spectrum, A^^(fe) and a smoothing window function, fk, 

cr^ = j Al^{k)\h\"^ . (CI) 

We use the recipe outlined in Peacock (1999) to derive 
the nonlinear power spectrum from the linear one. For the 
smoothing window we choose a Gaussian function /^^ = 
(Hui & Gnedin 1997, Nusser 2000), where xj is 
the Jeans length scale. 

The clumping factor can be defined as 

/ciu..p = ||| = {(l + 5)^> = l + {5^>, (C2) 

where p = (p)(l + (5) is the gas density at any point in the 
universe, (p) is the mean gas density, and 5 is the density 
fluctuation. 
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